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ABSTRACT 

A new method of imposing boundary conditions in the pseudospectral approximation of 
hyperbolic systems of equations is proposed. It is suggested to collocate the equations, not 
only at the inner grid points, but also at the boundary points and use the boundary conditions 
as penalty terms. In the pseudospectral Legrendre method with the new boundary treatment, 
a stability analysis for the case of a constant coefficient hyperbolic system is presented and 
error estimates are derived. 
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1. Introduction 


The importance of the correct numerical implementation of the boundary conditions, 
in approximating hyperbolic systems of equations, is widely recognized. The pioneering 
works of Gustafsson, Kreiss & Sundstrom [6] and Osher [7], provide the stability theory 
for the boundary conditions treatment, in the framework of the finite differences method. 
Basically, the conditions for stability are reduced to an algebraic problem. 

The role of boundary conditions in spectral or pseudospectral methods is even more 
crucial than in finite differences method. One of the reasons is that, since spectral methods 
are global methods, the behavior at one point affects the computation in the whole domain, 
so that the information at the boundary propagates very fast. Therefore, if on one hand the 
spectral algorithms do not necessitate special treatment at the boundary (thus in general 
there are no numerical boundary conditions ) , on the other hand a non correct specification 
may cause explosive instabilities. 

In this paper we introduce a new method of applying boundary conditions when 
approximating hyperbolic systems of equations by the collocation method. The novel idea 
is to collocate the differential system at all the grid points (included the boundary points) 
and use the boundary conditions as a penalty term. The same idea was developed in [3] 
for the pseudospectral Chebyshev discretization of a scalar hyperbolic equation. Here we 
prove stability and show convergence estimates for the pseudospectral method based on 
the Legendre nodes. 

The paper is structured as follows. In section 2 we describe the method in the scalar 
case and prove error estimates for the pseudospectral Legendre method. In section 3 we 
treat the system in the diagonal form with coupling through the boundary conditions. A 
complete convergence result is shown. In section 4 we show how to implement the method 
in a more general case. 
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2. Estimates for the scalar equation 

In this section we derive error estimates for the new pseudospectral Legendre approxima- 
tion to a scalar hyperbolic constant coefficients problem. Let U = U(x,t) be the solution 
to: 




'W 

V t = aU x 

x 6 [-1,1], 0 <t <T, 

(2.1) 

i 

(ft) 

17(1, t) = h(t) 

0 < t < T, 



<(c) 

U(x, 0 ) = f(x) 

*€[-1,1], 


where h and / are given functions and aGR, a > 0. 

Let N to be an integer. In the pseudospectral Legendre method we approximate U 
by v = t;(x,i), which is required to be a polynomial of degree at most N in the variable 
x for any 0 < t < T. This is done by demanding that t; satisfies equation (2.1) (a) at the 
grid points xy (j — 1, TV). The points xy (j = 0 ,N) are the nodes of the Gauss-Lobatto 
quadrature formula and the extrema in [-1,1] of the N th degree Legendre polynomial Pj <j. 
More precisely we take xo = 1, xn = —1, while xy (j = 1 ,N — 1) are the zeroes of P' N in 
decreasing order. 

The choice of this particular grid allows an accurate evaluation of integrals by summing 
over the grid values. Namely, let wy (j — 0, N) be the weights of the Gauss-Lobatto 
formula, then for any polynomial p of degree at most 2 N — 1 the following equality holds 
(see, for instance [2]): 


( 2 . 2 ) 


,1 N 

/ Pdx = ^2 P( x i)U} 

J ~ 1 3-0 


In the following, we will set u = uq = u> ^ • 

The new method involves a different treatment of the boundary condition (2.1) (6) (see 
[3]), rather than imposing exactly such a condition. In fact, we are concerned in finding v 
such that: 


(2.3) 


r (a) v t = av x 

< (6) Vt(l,t) = at>*(M) -7(v(l,t) - h(t)), 
,(c) V(xy,0) =/(xy) 


at x = Xy, j — l,N, 


j - 0, N. 


The coefficient 7 > 0 in (2.3) ( b ) will 


be specified later for the analysis of stability. 
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It is convenient to compare v(-,t) with some projection of the solution ?/(•,£) in the space 
of polynomials. To this end we introduce two projection operators denoted by In and IlAr 
respectively. 


Definition 2-1 InU is the polynomial of degree at most N that interpolates U at the 
points Xj (j = 0 ,N), i.e.: InU(xj) = U(xj), j =0,N. 

Let H a {- 1,1), o e R be the usual Sobolev spaces with if 0 (-1,1) = L 2 (-l,l). 

Definition 2-2 ITjvJ/ is the polynomial of degree at most N that is the best approximation 
of U in the if 1 (—1,1) norm with the condition: IIjvC 7’(±1) = I7(±l). 

The two following error estimates concerning the projectors In and 11^ can be found 
in [1]: 

(2.4) || U - InUWhu-ui) < CN^WUWh^-^), 

WU e H a {- 1, 1), °>\, 0 < /i < a; 

Zt 

(2.5) II u - n*lf||tfM(-i,i) < CN»-°\\U\\ H ', [ - UXh 

VU G 1), 0>1, 0 <n<0. 

In the next theorem we estimate the error 6 = II nU — v, where U is the solution to 
(2.1) and v is the solution to (2.3). 

Theorem 2-1 

Let 'i > — then we have: 

f — fjj 


( 2 . 6 ) 


y\6 2 (xj,T)uj + a f ($ 2 (l,t) + S 2 (-l,t))dt < 

T^o Jo 


< e 


N rT N \ 

T(I N f - n N f) 2 {xj)uj + T / , 

J = 0 J ° 3=0 J 


where Q — a[ntft/x — (ITmI/)*]. 
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Proof. Applying ITjv to (2.1) (a), one gets: 


(2.7) 


r (n^C/) t = a(n^Z7) a: + Q at x = x 3 , j = 0,N, 
< ( ITjv C/ ) ( 1 , t) = h(f) 0 < t < T, 

, (n*t/)(xy,o) = (n N f){x 3 ) j = o, n. 


Therefore, upon subtracting (2.3) from (2.7), we have the following error equation: 


( 2 . 8 ) 


(а) <5 t = a6 x + Q at x = x 3 , j =- 1 ,N, 

(б) 6 t (M) = «MM) + <5(1, t) - 7<5(l,<), 

(c) 5(xy,0) = (TIw/ — In/){xj) 3 = 0, N. 


Note that 8 is a polynomial of degree at most N . 

Multiplying (2.8)(a) by 6(x 3 ,t)u>j ( j = 1 ,7V), (2.8)(6) by 5(l,f)w 0 , and summing, we get 
(using (2.2)): 


N 

(2.9) 5^(5 t^)(*y,O w j = a 

y=o 


i: 


8 x 8dx + 


N 

+ ^(Q<5)(xy,f)wy - iw6 2 (l,t), 0 < t < T. 

3 = 0 

Integration by parts yields: 


N 


N 


(2.10) — ^26 2 (xj,t)uij = (a - 27 w)< 5 2 (M) -aS 2 (-l,t) + 2^{Q8)(x j ,t)u 3 < 


i=o 


3=0 
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1 1 II i 1 1| IIRIII II 


N N 

< -a( 6 2 (l,t) + 5 2 (-l,f)) + ™ + r SG a (*y»O w /» 0 ^ 1 ^ T > 


3=0 


3=0 


where we used the relation a — 2 ')u> < —a. Finally the Gronwall lemma yields (2.6).0 

Recall now that there exist two positive constants C i and C% such that one has (see 
[1], page 286): 


( 2 . 11 ) 


r 1 N f 1 

Ci / p 2 dx < \^p 2 (xj)ujj < C 2 / p 2 dx, 
J ~ l 3=0 J ~ l 


for any polynomial of degree at most N. So, we can prove the following result. 

Theorem 2-2 

Let U e L 2 {Q,T\H <T {-l,\))C\L 00 {0,T\H a - 1 {-\,l)),o > 1 be the solution of { 2.1). Let v 
be the solution of (2.3) with 7 > jj . Then we have the error bounds: 


( 2 . 12 ) 


|(u -»)(•, rjiu.,.,,!, < 


< 


CN 1 ” (yf ||H|jL s (0,Ti/f°(-l.l)) + l|P(-.r) III,— !(_,,!) + ||/||h.->(-1,1)) i 


(2.13) 




< CN 1 ° (Vf ||t/||z,3(0,T;ir‘ J (-l,l)) + l|/||/f'’- 1 (-lil)) • 

Proof. Note first that by (2.5) Q can be estimated in the following way: 

(2.14) ||Q||l*(-i,i) < + a\\U N U x — U x ||l j (-i,i) < 


< CN l ~ a 




a > 1. 
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Then, (2.12) is easily obtained by the previous theorem, by (2.11), (2.4) and by the 
triangle inequality: 


(2.15) \\U - t/||i,a(_i,i) < II U - ITivZ7||x,a(— 1,1) + \\n N U - v||x,a ( — j.,1) • 


To show (2.13) it is sufficient to observe that, by definition, JInU coincides with U at 
x = ±1 for any t , so that: 


(2.16) 



U - U) 2 {±l,t)dt 


= 0 , 


and then use the estimates shown above. O 


Note that the same results hold when a < 0 and the role of the boundary points x = 1 
and x = — 1 is interchanged. We note also that, as in the Chebyshev case (see [3]), 7 has 
to be proportional to N 2 in order to prove a stability result. 
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3. Error analysis for a system of equations 

Let ki,k 2 be two integers and denote by A ^ the following ki x ki (t = 1,2) diagonal 
matrices: 


A W = diag{a^\ • • 

where a >0, j = l,fci and ay ^ > 0, j = 1,&2- 

Let R be a k 2 x ki matrix and L a k% X &2 matrix. Let r and / be the norms of the 
operators R and L respectively, i.e.: 

r = ll#ll£(R fc i,R''J)> 1 ~ ll^ll£(R fc >,R fc i)- 

In the following we assume that r and l satisfy the relation: 


(3.1) 


0 < rl < 1. 


Denote by < ■ , • >» and by || • ||i the scalar product and the norm in R fc ’ (t = 1,2) 
r espect ively. For a given positive definite ki x ki diagonal matrix M ( ‘), we shall denote by 
V9) the corresponding positive definite fcj x ki diagonal matrix whose entries are the 
square roots of the entries of . Note that we have the equalities. 


(3.2) < >i = < >i - < >* > 


* = i,2. 

We are now ready to state the differential problem. ^Ve are concerned with finding 
the functions vector U = U(x,t ), where U = and U^(x,t) G R Vr G 

[-1,1], Vt G [0,r], i = 1,2; such that: 


(3.3) 


r £/ t (1) = -A^U { x l) 

< u\ 2) = AWui 2 \ x G [—1, 1], 0 < t < T, 
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(3.4) 


trt 1 >(-i,0 = £tri , )(-i,0 + »iW 

C/< 2 >(1,<) =RUW(l,t) + g 2 (t), 0 < t < T, 


(3.5) 


ir(‘>(x,0) = f l (x) 

UW(x,0) = f 2 (x), x £ [—1,1], 


where fi,gi (t = 1,2) are given functions. 

In the Legendre collocation method we seek u = u(x,t ) with u = (u^,^ 2 )), u^(x,t) G 
R fcf , Vx G [—1,1], Vt G [0,r], t = 1,2. The vectors i = 1,2, whose components 
are polynomials in the variable x of degree at most N, are determined according to the 
following collocation scheme: 


(3.6) 


I u t (1) = at x = xy, j = 0,N - 1, 

| = A^ui 2 ^ at x = xy, j — 1 ,N, 


u t — 1*0 — — (— 1, t) — r( 1 ){u( 1 )(-l,t) — Lu^(— 1, t) — ^i(t)}, 

«{ 2) (1»0 = ^ 2) ui 2) (l,<) - T( 2 ){u( 2 )(l,t) - RuW(l,t) - g 2 (t)}, 


('u( 1 )(x i , 0) = / 1 (x i ) j = 0,N, 
(3.8) l 

(u( 2 )(xy,0) = / 2 (xy) i = 0,J\T. 


Here TW, * = 1,2 are positive definite ki x k{ diagonal matrices to be specified later. 
Note that we took into account at the points x 0 = 1 and x# = 1 both the equations and 
the boundary conditions. 
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In order to study the convergence of u to U when N tends to +oo, we introduce another 
collocation scheme. Namely, we look for v = v(x,t), v = (u (1 J ), such that: 


(3.9) 


at x = Xj, 
v t (2) = A^v i 2) at x = Xj, 


j = 0, JV — 1, 

j = hN, 


„W(-i ,o = -A^ti^-M) - rU){ v (i)(_i,t) - 

„t a >(M) = A^v ( x 2 \l,t)-T^{v^(l,t)-U^(l,t)}, 


(vW(xj,0)=h(xj) j = 0,N, 

(3.11) { 

{vW(xj,0) = Mxj) j = 0,N. 

This time the system is totally uncoupled. Instead of comparing U and u we estimate 
the error between u and v. The scalar analysis will enable us to estimate the difference 

U- v. 

Theorem 3-1 

Let then we have the estimate: 


(3.12) (r ||2^ 1 W 1 )(xj ) T , )|| 2 + / \\B^w^(xj,T)\\^j Wj < 

3=0 

< — r/ ^_ y J* | H(t/ {1) - V (1) )(l,t)||f + \Jl ||(tf (2) - » (2) )(-l,0llaj dt > 

where w — u — v and * , i = 1*2. 

Proof. From (3.6)-(3.11) it is clear that w satisfies the equations: 
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at x = Xj, j = 0,N - 1, 
atz = x y , j = 1,7V, 

(3.14) 

«; t (1) (-l,t) = -^ 1 )4 1) (-l,0-r( 1 ){ U ;( 1 )(-l,t)-X«,(2)(-i ) t) + L £ (2)(_ ljf )j } 
w t (2) (l,0 = A^wi 2) (l,t) -T( 2 ){«;( 2 )(l,t) - JZwWfl.f) +ReW(l,t)}, 


j Wf 1 ^ = —A^wi 1 ^ 
1 w { t 2) = AWwl 2) 


(3.15) 


^(xy.O) =0 j = 0,N, 

u;( 2 )(xy,0) =0 j = 0, N, 


where t/*), t — 1,2. Note that (3.4) has been used to eliminate U^(— 1,<) 

and UW{l,t). K ’ } 

Now, let £>(*), t = 1,2 two diagonal positive definite x matrices to be specified 
later. Multiplying the first set of equations in (3.13) and (3.14) by {D^w^){xj,t)uj and 
summing up on j = 0, N, we get by (2.2): 


N x 

(3.16) Y, < >, (*y,«)u>, = - / < aMwW.dWwW >, (x,t)dx + 

y=o j - 1 


-w<r( 1 ) U ;( l )-r{ 1 )W 2 ), J D( 1 )«;( 1 ) >! {-l,t)-u<T^Le^,D^w^ > x (-1 ,<) < 


< >! (1 ,t) 


+ ^ >1 (— l,f) + 
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-w < I^V 1 ) - r (1 )W 2) ,P (1) u; (1 > >1 (-1,4) + 


+ <tMwW,dMwM >1 (-1,4) + ££ (2) (— 1,4)||?, 

2 .2t] 

where 0 < r/ < 2. In the same way we have: 


N 


(3.17) 


D < >2 (*/,4)wy < 


i=o 



< a( 2 W 2 ),£>( 2 W 2 ) >2 (1,4) - ^ < a( 2 W 2 ),Z?( 2 W 2 ) >2 (—1,4) + 

u 


u < TWwW - tWRw( 1 \dWwW > 2 (I,t) + 


+ < r( 2 V 2 \D (2 V 2 > > 2 ( 1 , 4 ) + ^-||>/i?( 2 )r( 2 ) JRe (1) (i,4)|g. 

2 217 

Note that 77, D (*) and Z)( 2 ) are still to be specified. We start by setting r} = 2(1 — Vrl). 
Recalling the hypotheses on I'M, by summing (3.16) and (3.17) one obtains: 


(3.18) 


^2 [||\^W«; (1) ||? + ||\/.D( 2 )u; ( 2) || 2 ] (zy,4)wy = 


J=o 




-l=<AWRwl l \DWwW > 2 +||\/a( 2 )P( 2 )u;( 2 )|| 2 1(1,4) + 
vW J 
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-[||\/a(‘)£)(‘)w< 1 »||;--^= < aWLwW,dMw( 1 > >! +||\/A( 2 )D( 2 )t»( 2 >|| 2 ](-l,i) + 


+ 2(1 _^ )v ^ (IIv'dWAW J?£ (1) (l,t)||| + Hv'cCUU) ie< 2 >(-l,t)ll?) • 

At this point we want the two first terms in brackets on the RHS of (3.18) to be 
positive. This is true if we choose: 

0(‘)=(Af I »)'‘r 2 , D< 2 > = (a< 2 >)'\/. 

In fact, the first of those terms becomes: 

r 2 ||u/W||? — 2\/ri < Rw^\w^ >2 H- r/||tt;f 2 ^ |||, 

and this is positive since: 


\2y/rl < RwW,wW > 2 I < \\RwW\\l + r/||^ 2 )|| 2 <r 2 ||u>( 1 )|| 2 + W||u;( 2 )|| 2 . 

Similar arguments hold for the second term. This allows us to write the following 
inequality: 


(3.19) 


d_ 

dt 


y=o 


+ 





< 


- ( < n jhtI 1 ( 1 -« ) ni + ^(-mik) < 
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lit HI INI 


* (Vt ||£< 1,(m)|i? + V> (2,( - M)iii ) ■ 

Finally, by integrating in time, we obtain (3.12).0 

Using the results of the previous section we can finally prove our main convergence 
theorem. 

Theorem 3-2 

Let U be the solution of (3.3)-(3.4)-(3.5) and let u be the solution to (3.6)-(3.7)-(3.8) with 
r<‘> = 4^ • Then one has: 


(3.20) 


E n ^ 0 



< 

L*(-l,l) 


< CN l ~ a T 



? 

C°([0,T];J/ <, ( — 1,1)) 


where C only depends on l,r,cr and 4^), * — 1,2. 

Proof. We first write: U — u — U — v — w. Then a bound for each component of 
£/U) — t;W, * = 1,2 is given by (2.12), while for w we have: 


(3.21) 


X> (i) lin (*,r) 


1 


L’l- hi) 


^||[pW] V^ll?) (;T) 


< 


\i = 1 


IX.*(-1,1) 
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< max 
*= 1,2 


£(R*i,R*.) 


5^||BWwW||? (;T) 




< C[r,l) max {(aj.^) 2 } ||^ (1) -« (1) ||?(1 ,t)dt + ||*/ (2) - v (2) \\l(-l,t)dt 


where we used theorem 3-1. Finally the last term in (3.21) is estimated as in (2.13).0 
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4. Suggestions for the implementation of non diagonal systems 

In this section we discuss the implementation of our new approach to the general hyperbolic 
system: 


(4.1) 


Zt — H Z x , 


where Z = Z(x,t) is a k component vector and H is a constant coefficients k x k matrix 
with ki negative eigenvalues and k 2 positive eigenvalues (k = ky + k 2 ). 

The following boundary conditions are imposed at x = —1: 

(4-2) (*‘ 

where B n is a ^ x k t matrix, B l2 is a fci x k 2 matrix and h x is a given components 

vector. Besides, the following boundary conditions are imposed at x = 1: 


(4.3) 



bL) 2f(M) 


where B 2i is a k 2 x fcj. matrix, £22 is a fc 2 x k 2 matrix and h 2 is a given k 2 components 
vector. 

Suppose that there exists a nonsingular matrix: 



where T i3 is a k { x k 3 matrix, such that the change of variables Z - TU diagonalizes the 
system (4.1). Thus, we get: 


(4.4) 


dU 

dt 


T~ l HT 


dU 

dx 


(-A™ o \d£ 
^ 0 AW ) dx ’ 


where U and A^ % \ 


t = 1,2 have been defined in section 3. 
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The boundary conditions are respectively transformed as follows: 


(4.5) (BnTn + B l 2 T 2 l )UM{-l,t) + {B tl T 12 + B l 2 T 22 )U^{-l,t) = h t {t), 

(4.6) {B 2 iTu + B 22 T 2 i)U^^ (1, t) + ( B 2 \Ti 2 + B 22 T 22 )U^(l,t) = h 2 (t). 


Therefore (4.2) and (4.3) are equivalent to (3.4) and (3.5) if and only if the matrices 
+ B 12 T 21 and B 2 iT i2 + B 22 T 22 are invertible. In this case we have: 

L = — (BuTn + B12I21) 1 (Bn 2 ’i 2 + B12T22), 


R — — (B21T12 + B 22 T 22 ) 1 (B 2 iTu + B 22 T 2 i), 


9 i — (BiiTii + B12T21) 1 h i , 


92 = (B21T12 + B 22 T 22 ) 1 h 2 . 

We would like to show how to apply the scheme (3.6)-(3.8) directly to the system 
(4.1)-(4.3). Assuming the hypotheses of theorem 3-2 we have T^*) = /SA^ l \ i = 1,2, where 
= N(N + l)/2y/ri. We also set (denoting by u the approximation to U ): 

/ « (1) (z> t) - W 2 ) (z, t ) - g t ( t ) \ 

E(x,t) = ( 1 = 

V « (2) (z, 0 - Bu W (z, t) - g 2 (<) ) 

= (-R 

Thus (3.6)-(3.7) can be written in the form: 


(4.7) -^(x,-,t) = A-^(x 3 -,t) + P 6 N 3 A^E{-I,t) - p 6 0 j A^E{l,t), j = 0,N, 
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where ^ , A<+> = ^( 2 ) ) » and A = A* ) + A(+>. 

Let us now define z = Tu. The function z satisfies: 


(4.8) | £(*,,«) = ff + fiS Hi (TA^T- l )TB(-l,t)+ 

-pSojiTA^T-^TZil^), j = 0,iV. 

Defining B = T ^ j y T -1 , we get: 

(4.9) 

Finally, taking fl**) = TA^T -1 and substituting in (4.8), we obtain the pseudospectral 
scheme to approximate (4.1)-(4.3), namely: 

(4.10) ~(T„t) = nf x (x„t) + P6 Ni H'-> [b*(-1.0 - r (£) w] + 

-/Mo ,H< + > B*(l,t)-r^)(t)], 1=0, J V. 

This is equivalent to collocate the equation (4.1) at all the points with some suitable 
penalty terms, deriving from the boundary conditions, added at the points x = ±1. It is 
clear that the same convergence estimates of theorem 3-2 also apply for the error Z -z = 
= T{U-u). 
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